## VAA effects meta-analysis ----------------------
## meta-analysis ----------------------------------

# load packages, data, and functions
source("packages.R")
load("vaa_df_prep.RData")
source("functions.R")

#covert NA to ?
vaa_turnout$vaa_usage_sample <- gsub("NA", "?", vaa_turnout$vaa_usage_sample)
vaa_vote$vaa_usage_sample <- gsub("NA", "?", vaa_vote$vaa_usage_sample)
vaa_knowledge$vaa_usage_sample <- gsub("NA", "?", vaa_knowledge$vaa_usage_sample)

# count number of effects
sum(
  nrow(vaa_turnout),
  nrow(vaa_vote),
  nrow(vaa_knowledge)
)

############################################################################
## Meta-analysis -----------------------------------------------------------
############################################################################

## random-effects models  -------------

# political knowledge
knowledge_raneff <- rma(yi = pcor, 
                        sei = pcor_se, 
                        data = vaa_knowledge,
                        slab = study_label,
                        method = "REML")
summary(knowledge_raneff)

## multilevel meta-analysis with cross-classified random effects -------------

# turnout
turnout_mleff <- rma.mv(yi = log_odds, 
                        V = log_odds_se^2, # square standard errors to feed sampling variances
                        random = list(~ 1 | study_label, ~ 1 | election),
                        data = vaa_turnout,
                        slab = study_label,
                        method = "REML",
                        tdist = TRUE)
summary(turnout_mleff)
compute_I2(data = vaa_turnout, model = turnout_mleff, vi = "log_odds_se")



# vote
vote_mleff <- rma.mv(yi = log_odds, 
                     V = log_odds_se^2, # square standard errors to feed sampling variances
                     random = list(~ 1 | study_label, ~ 1 | election),
                     data = vaa_vote,
                     slab = study_label,
                     method = "REML")
summary(vote_mleff)
compute_I2(data = vaa_vote, model = vote_mleff, vi = "log_odds_se")

## Figures 3-5. Forest plots with effect sized (odds ratios) of VAA usage

##forest plot with study type
vaa_turnout <- vaa_turnout %>% mutate(design_short = case_when(vaa_turnout$study_design == "Experimental" ~ "EXP",
                                                               vaa_turnout$study_design == "Observational - no panel" ~ "O-NP",
                                                               vaa_turnout$study_design == "Observational - panel" ~ "O-P",
                                                               vaa_turnout$study_design == "Observational - selection models" ~ "O-SM"
))

vaa_vote <- vaa_vote %>% mutate(design_short = case_when(vaa_vote$study_design == "Experimental" ~ "EXP",
                                                         vaa_vote$study_design == "Observational - no panel" ~ "O-NP",
                                                         vaa_vote$study_design == "Observational - panel" ~ "O-P",
                                                         vaa_vote$study_design == "Observational - selection models" ~ "O-SM"
))

vaa_knowledge <- vaa_knowledge %>% mutate(design_short = case_when(vaa_knowledge$study_design == "Experimental" ~ "EXP",
                                                                   vaa_knowledge$study_design == "Observational - no panel" ~ "O-NP",
                                                                   vaa_knowledge$study_design == "Observational - panel" ~ "O-P",
                                                                   vaa_knowledge$study_design == "Observational - selection models" ~ "O-SM"
))


pdf(file = "forest-turnout-ml-type.pdf", height = 7, width = 8, family = "URWTimes")
par(oma = c(1,0,0,0))
par(mar = c(1.5,3,0,.5))
header_line <- nrow(vaa_turnout) + 2
forest.rma.mod(turnout_mleff,
               slab = vaa_turnout$study_label,
               xlim = c(-14.5, 6), 
               at = log(c(0.25, 1, 5, 12)), 
               atransf = exp, 
               showweights = TRUE,
               ilab = cbind(vaa_turnout$design_short ,vaa_turnout$country, vaa_turnout$election_year, vaa_turnout$sample_size, vaa_turnout$vaa_usage_sample), 
               ilab.xpos = c(-10.3, -6.5, -5, -3.5, -2), 
               ilab.pos = c(4, 2, 2, 2, 2),
               cex = 0.75,
               xlab = "",
               mlab = "Overall estimate (CCREM with study and election REs)")
op <- par(cex = 0.75, font = 2) 
text(-14.5, header_line-.2, "Author(s) and Year", pos = 4) 
text(c(-10.3, -6.5, -5, -3.5, -1.8, 1.5, 3.8, 6),  header_line-.2, c("Type"," Country","Year", "Size", "VAA users", "Observed outcome", "Weight", "OR [95% CI]"), pos = c(4, 2, 2, 2, 2, 2, 2, 2))
text(c(-6.5, -3.25, 2.25), header_line + 2.5, c("Election", "Sample", "Estimates"), pos = c(1, 1))
arrows(x0 = -8, x1 = -5.25, y0 = header_line + .75, length = 0)
arrows(x0 = -4.5, x1 = -2, y0 = header_line + .75, length = 0)
arrows(x0 = -1.5, x1 = 5.75, y0 = header_line + .75, length = 0)
par(op)
dev.off()

pdf(file = "forest-vote-ml-type.pdf", height = 5, width = 8, family = "URWTimes")
par(oma = c(1,0,.5,0))
par(mar = c(1.5,3,.75,.5))
header_line <- nrow(vaa_vote) + 2
forest.rma.mod(vote_mleff,
               slab = vaa_vote$study_label,
               xlim = c(-14, 6), 
               at = log(c(0.25, 1, 5, 12)), 
               atransf = exp, 
               showweights = TRUE,
               ilab = cbind(vaa_vote$design_short, vaa_vote$country, vaa_vote$election_year, vaa_vote$sample_size, vaa_vote$vaa_usage_sample), 
               ilab.xpos = c(-10, -6.5, -5, -3.5, -2), 
               ilab.pos = c(4, 2, 2, 2, 2),
               cex = 0.75,
               xlab = "",
               mlab = "Overall estimate (CCREM with study and election REs)")
op <- par(cex = 0.75, font = 2) 
text(-14, header_line-.2, "Author(s) and Year", pos = 4) 
text(c(-10, -6.5, -5, -3.5, -1.8, 1.5, 3.8, 6),  header_line-.2, c("Type","Country", "Year", "Size", "VAA users", "Observed outcome", "Weight", "OR [95% CI]"), pos = c(4, 2, 2, 2, 2, 2, 2, 2))
text(c(-6.5, -3.25, 2.25), header_line + 2.5, c("Election", "Sample", "Estimates"), pos = c(1, 1))
arrows(x0 = -8, x1 = -5.25, y0 = header_line + .75, length = 0)
arrows(x0 = -4.25, x1 = -2, y0 = header_line + .75, length = 0)
arrows(x0 = -1.25, x1 = 5.75, y0 = header_line + .75, length = 0)
par(op)
dev.off()

pdf(file = "forest-knowledge-type.pdf", height = 2, width = 8, family = "URWTimes")
par(oma = c(1,0,0,0))
par(mar = c(1.5,3,1,.5))
header_line <- nrow(vaa_knowledge) + 2
forest.rma.mod(knowledge_raneff,
               slab = vaa_knowledge$study_label,
               xlim = c(-1.2, .68), 
               at = c(-.1, 0, .1, .2, .3), 
               ilab = cbind( vaa_knowledge$design_short, vaa_knowledge$country, vaa_knowledge$election_year, vaa_knowledge$sample_size, vaa_knowledge$vaa_usage_sample), 
               ilab.xpos = c(-0.9,-.55, -.43, -.3, -.13), 
               ilab.pos = c(4, 2, 2, 2, 2),
               cex = 0.75,
               xlab = "",
               showweights = TRUE,
               digits = 2,
               psize.shrinkage = .3,
               psize.baseline = 1.25,
               mlab = "Random-effects model",
               top = 4)
op <- par(cex = 0.75, font = 2) 
text(-1.2, header_line-.2, "Author(s) and Year", pos = 4) 
text(c(-.9, -.55, -.43, -.3, -.13, .22, .45, .68),  header_line-.2, c("Type", "Country", "Year", "Size", "VAA users", "Observed outcome", "Weight", "P. corr. [95% CI]"), pos = c(4, 2, 2, 2, 2, 2, 2, 2))
text(c(-.58, -.25, .3), header_line + 2.25, c("Election", "Sample", "Estimates"), pos = c(1, 1))
arrows(x0 = -.7, x1 = -.46, y0 = header_line + .75, length = 0)
arrows(x0 = -.37, x1 = -.15, y0 = header_line + .75, length = 0)
arrows(x0 = - .03, x1 = .66, y0 = header_line + .75, length = 0)
par(op)
dev.off()


# show where weights come from --------------------

# extracted weights
weights(turnout_mleff) 

# extracted weights scaled in percent
wi <- diag(chol2inv(chol(turnout_mleff$M)))
wi100 <- wi/sum(wi) * 100
plot(wi100, weights(turnout_mleff))
abline(0,1)

